Code to run analysis and generate figures for McNeall et al. (2023) “Constraining the carbon cycle in JULES-ES-1.0”
source("JULES-ES-1p0-common-packages.R")
Loading required package: spam
Loading required package: dotCall64
Loading required package: grid
Spam version 2.7-0 (2021-06-25) is loaded.
Type 'help( Spam)' or 'demo( spam)' for a short introduction
and overview of this package.
Help for individual functions is also obtained by adding the
suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
Attaching package: ‘spam’
The following objects are masked from ‘package:base’:
backsolve, forwardsolve
Loading required package: viridis
Loading required package: viridisLite
See https://github.com/NCAR/Fields for
an extensive vignette, other supplements and source code
Loading required package: lhs
Loading required package: parallel
Loading required package: viztools
source("JULES-ES-1p0-common-functions.R")
source("JULES-ES-1p0-common-data.R")
getStandardMemberLongName <- function(ensloc, variable){
# Helper function to get the long name of output variables (and the units)
ensmember <- 'S3'
fn <- paste0(ensloc,'JULES-ES-1p0_',ensmember,'_Annual_global.nc')
nc <- nc_open(paste0(fn))
ln <- nc$var[[variable]]$longname
un <- nc$var[[variable]]$units
return(list(longname = ln, units = un))
}
# for(i in y_names_select){
#
# out <- getStandardMemberLongName(ensloc_wave00, i)
# print(i)
# print(out)
# }
low_npp_ix <- which(Y[,'npp_nlim_lnd_sum'] < 1e5)
# code from https://stackoverflow.com/questions/28182872/how-to-use-different-sets-of-data-in-lower-and-upper-panel-of-pairs-function-in
#X <- matrix(runif(300), ncol=3)
#Y <- matrix(c(sort(runif(100, 0, 10)),
# sort(runif(100, 0, 10)),
# sort(runif(100, 0, 10))), ncol=3)
#pdf(file = 'figs/fig02.pdf', width = 12, height = 10)
#pdf(file = 'figs/run-failure-pairs.pdf', width = 12, height = 10)
x1 <- X[low_npp_ix, ]
x2 <- X_nlevel0
XY <- rbind(x1, x2)
pairs(XY,
lower.panel=function(x, y, ...) {
Xx <- x[seq_len(nrow(x1))] # corresponds to X subset
Xy <- y[seq_len(nrow(x1))] # corresponds to X subset
#usr <- par("usr"); on.exit(par(usr))
#par(usr = c(range(x1[, -ncol(x1)]), range(x1[, -1]))) # set up limits
points(Xx, Xy, col = zblue, pch = 19, cex = 0.8)
# if(par('mfg')[2] == 1) axis(2) # if left plot, add left axis
#if(par('mfg')[1] == ncol(x1)) axis(1) # if bottom plot add bottom axis
},
upper.panel=function(x, y, ...) {
Yx <- x[(nrow(x1) + 1):length(x)] # Y subset
Yy <- y[(nrow(x1) + 1):length(y)] # Y subset
#cntr <- outer(Yx, Yx, FUN='*') # arbitrary function for contour
# usr <- par("usr"); on.exit(par(usr))
#par(usr = c(range(x2[, -1]), range(x2[, -ncol(x2)]))) # set up limits
points(Yx, Yy, col = zred, pch = 19, cex = 0.8)
#contour(Yx, Yy, cntr, add=TRUE)
#if(par('mfg')[2] == ncol(x2)) axis(4) # if right plot, add right axis
#if(par('mfg')[1] == 1) axis(3) # if top plot, add top axis
},
#tick=FALSE, # suppress the default tick marks
#line=1,
gap = 0,
xlim = c(0,1), ylim = c(0,1),
labels = 1:d,
oma = c(2, 18, 2, 2)) # move the default tick labels off the plot
reset()
legend('left', legend = paste(1:d, colnames(lhs)), cex = 1.1, bty = 'n')
legend('topleft', pch = 19, col = c( zred, zblue), legend = c('failed', 'zero carbon cycle'), bty = 'n', inset = 0.02, cex = 1.1 )
#dev.off()
Are there hard thresholds after which the model always fails? (Question from review comments). Run failures are in t
y_level0 <- Y_level0[,'npp_nlim_lnd_sum']
pdf(file = 'figs/marginal_failures.pdf', width = 10, height = 10)
par(mfrow = c(5,7), mar = c(3,1,3,1), oma = c(1,1,5,1))
for(i in 1:p){
plot(X_level0[,i], y_level0, xlab = '', ylab = '', main = colnames(X_level0)[i], xlim = c(0,1))
points(X_nlevel0[,i], rep(0, nrow(X_nlevel0)), col = zred, pch = 19)
points(x1[, i ], rep(0,nrow(x1)), col = zblue, pch = 19)
}
mtext('NPP', outer = TRUE, side = 3, cex = 2, line = 2)
reset()
legend('bottomright', pch = 19, col = c('black', zred, zblue), legend = c('ran','failed', 'zero carbon cycle'), bty = 'n', inset = 0.05, cex = 1.1 )
dev.off()
null device
1
## Is there a pattern to the wave01 timeseries outliers?
pdf(file = 'figs/X_wave01_timeseries_outliers.pdf', width = 12, height = 12)
pairs(X_wave01_train[ts_outliers_ix_wave01, ],
gap = 0,
xlim = c(0,1),
ylim = c(0,1)
)
dev.off()
null device
1
It’s important to remember that the design of the experiment is multiplication factors of the original parameters. This might be important for the “hold” value in a sensitivity analysis, as the “standard” value and the median value of the ensemble will not be the same.
#pdf(file = 'figs/lhs_range.pdf', width = 6, height = 8)
#pdf(file = 'figs/fig01.pdf', width = 6, height = 8)
par(las = 1, mar = c(5,8,2,1))
lhs_min <- apply(lhs_wave0_wave01_all, 2, min)
lhs_max <- apply(lhs_wave0_wave01_all,2, max)
plot(lhs_max, 1:d, type = 'n', xlim = c(0,10), axes = FALSE, xlab = 'multiplying factor', ylab = '')
abline(v = 0, lty = 'dashed', col = 'grey')
abline(v = 1, lty = 'dashed', col = 'tomato2')
segments(x0 = lhs_min, y0 = 1:d, x1 = lhs_max, y1 = 1:d )
points(lhs_min, 1:d, pch = 20)
points(lhs_max, 1:d, pch = 20)
axis(2, at = 1:d, labels = colnames(lhs))
axis(1)
#dev.off()
Global mean for the 20 years at the end of the 20th Century. There is still a significant low bias on cVeg output.
wave00col <- 'skyblue2'
wave01col <- 'tomato2'
wave00col <- 'dodgerblue2'
wave01col <- 'firebrick'
rangecol <- 'grey'
# Histogram of level 1 constraints
hcol = 'darkgrey'
lcol = 'black'
#pdf(file = 'figs/fig05.pdf', width = 8, height = 8)
#pdf(file = 'figs/level_2_constraints_hists.pdf', width = 8, height = 8)
par(mfrow = c(2,2), fg = 'darkgrey', las = 1, oma = c(0.1, 0.1, 4, 0.1))
trunc <- function(x, vec){
dat <- x[x < max(vec) & x > min(vec) ]
dat
}
h <- hist(Y_const_level1a_scaled[,'nbp_lnd_sum'], main = 'NBP', xlab = 'GtC/year', col = makeTransparent(wave00col,150))
hist(trunc(Y_const_wave01_scaled [,'nbp_lnd_sum'], h$breaks) ,
col = makeTransparent(wave01col,150) , breaks = h$breaks, add = TRUE)
rug(Y_const_stan_scaled['nbp_lnd_sum'], lwd = 2)
polygon(x = c(0, 100, 100, 0), y = c(0, 0, 1000, 1000),
col = makeTransparent(rangecol, 60),
border = makeTransparent(rangecol))
h <- hist(Y_const_level1a_scaled[,'npp_nlim_lnd_sum'],col = makeTransparent(wave00col,150), main = 'NPP', xlab = 'GtC/year')
hist(trunc(Y_const_wave01_scaled [,'npp_nlim_lnd_sum'], h$breaks) ,
col = makeTransparent(wave01col) , breaks = h$breaks, add = TRUE)
rug(Y_const_stan_scaled['npp_nlim_lnd_sum'], lwd = 2)
polygon(x = c(35, 80, 80, 35), y = c(0, 0, 1000, 1000),
col = makeTransparent(rangecol, 60),
border = makeTransparent(rangecol))
h <- hist(Y_const_level1a_scaled[,'cSoil_lnd_sum'], col = makeTransparent(wave00col,150), main = 'Soil Carbon', xlab = 'GtC')
hist(trunc(Y_const_wave01_scaled [,'cSoil_lnd_sum'], h$breaks) ,
col = makeTransparent(wave01col,150) , breaks = h$breaks, add = TRUE)
rug(Y_const_stan_scaled['cSoil_lnd_sum'], lwd = 2)
polygon(x = c(750, 3000, 3000, 750), y = c(0, 0, 1000, 1000),
col = makeTransparent(rangecol, 60),
border = makeTransparent(rangecol))
h <- hist(Y_const_level1a_scaled[,'cVeg_lnd_sum'], col = makeTransparent(wave00col,150), main = 'Vegetation Carbon', xlab = 'GtC')
hist(trunc(Y_const_wave01_scaled [,'cVeg_lnd_sum'], h$breaks) ,
col = makeTransparent(wave01col,150) , breaks = h$breaks, add = TRUE)
rug(Y_const_stan_scaled['cVeg_lnd_sum'], lwd = 2)
polygon(x = c(300, 800, 800, 300), y = c(0, 0, 1000, 1000),
col = makeTransparent(rangecol, 60),
border = makeTransparent(rangecol))
reset()
legend('top', horiz = TRUE, fill = c(makeTransparent(wave00col, 150), makeTransparent(wave01col, 150), makeTransparent(rangecol, 60)), legend = c('Wave00', 'Wave01', 'AW range'))
#dev.off()
Just under a third. Points at a significant model discrepency in cVeg
Of the 400 members of the wave01 ensemble, 128 pass Andy Wiltshire’s Level 2 constraints.
length(level2_ix_wave01)
[1] 128
length(level2_ix_wave01) / ntrain_wave01
[1] 0.32
lcol_wave0 <- makeTransparent('dodgerblue2', 120)
lcol_wave01 <- makeTransparent('firebrick', 120)
lcol_wave01_level2 <- 'gold'
stancol = 'black'
linePlotMultiEns <- function(years, ens1, ens2, ens3, col1, col2, col3, ylab, main, ylim = NULL, ...){
# Plot wave00 and wave01 timeseries on top of one another
nt <- length(years)
if(is.null(ylim)){
ylim = range(c(ens1[,1], ens1[,nt], ens2[,1], ens2[ ,nt], ens3[,1], ens3[, nt]))
}
else ylim <- ylim
matplot(years, t(ens1), type = 'l', lty = 'solid',ylim = ylim, col = col1,
ylab = ylab, main = main, xlab = '',
bty = 'n', ...)
matlines(years, t(ens2), col = col2, lty = 'solid')
matlines(years, t(ens3), col = col3, lty = 'solid')
}
pdf(file = 'figs/fig03.pdf', width = 10, height = 12)
#pdf(file = 'figs/carbon-cycle-timeseries-waves-constrained.pdf', width = 10, height = 12)
par(mfrow= c(3,5), las = 1, mar = c(4,4,1,0))
linePlotMultiEns(years = years, ens1 = npp_ens_wave00[without_outliers_ix_wave00,],
ens2 = npp_ens_wave01[without_outliers_ix_wave01,],
ens3 = npp_ens_wave01[level2a_ix_wave01, ],
col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
ylab = 'GtC', main = 'NPP')
lines(years,npp_stan, col = stancol, lty = 'solid', lwd = 2)
linePlotMultiEns(years = years, ens1 = nbp_ens_wave00[without_outliers_ix_wave00,],
ens2 = nbp_ens_wave01[without_outliers_ix_wave01,],
ens3 = nbp_ens_wave01[level2a_ix_wave01, ],
col1 = lcol_wave0, col2 = lcol_wave01,col3 = lcol_wave01_level2,
ylab = 'GtC', main = 'NBP', ylim = c(-10,10))
lines(years, nbp_stan, col = stancol, lty = 'solid', lwd = 2)
linePlotMultiEns(years = years, ens1 = cSoil_ens_wave00[without_outliers_ix_wave00,],
ens2 = cSoil_ens_wave01[without_outliers_ix_wave01,],
ens3 = cSoil_ens_wave01[level2a_ix_wave01, ],
col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
ylab = 'GtC', main = 'cSoil', ylim = range(c(cSoil_ens_wave00[,1], cSoil_ens_wave00[,164])))
lines(years, cSoil_stan, col = stancol, lty = 'solid', lwd = 2)
linePlotMultiEns(years = years, ens1 = cVeg_ens_wave00[without_outliers_ix_wave00,],
ens2 = cVeg_ens_wave01[without_outliers_ix_wave01,],
ens3 = cVeg_ens_wave01[level2a_ix_wave01, ],
col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
ylab = 'GtC', main = 'cVeg')
lines(years, cVeg_stan, col = stancol, lty = 'solid', lwd = 2)
linePlotMultiEns(years = years, ens1 = lai_lnd_mean_ens_wave00[without_outliers_ix_wave00,],
ens2 = lai_lnd_mean_ens_wave01[without_outliers_ix_wave01,],
ens3 = lai_lnd_mean_ens_wave01[level2a_ix_wave01, ],
col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
ylab = 'GtC', main = 'Lai')
lines(years, lai_lnd_mean_stan, col = stancol, lty = 'solid', lwd = 2)
linePlotMultiEns(years = years, ens1 = rh_lnd_sum_ens_wave00[without_outliers_ix_wave00,],
ens2 = rh_lnd_sum_ens_wave01[without_outliers_ix_wave01,],
ens3 = rh_lnd_sum_ens_wave01[level2a_ix_wave01, ],
col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
ylab = 'GtC', main = 'RH')
lines(years, rh_lnd_sum_stan, col = stancol, lty = 'solid', lwd = 2)
linePlotMultiEns(years = years, ens1 = fLuc_lnd_sum_ens_wave00[without_outliers_ix_wave00,],
ens2 = fLuc_lnd_sum_ens_wave01[without_outliers_ix_wave01,],
ens3 = fLuc_lnd_sum_ens_wave01[level2a_ix_wave01, ],
col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
ylab = 'GtC', main = 'fLuc')
lines(years, fLuc_lnd_sum_stan, col = stancol, lty = 'solid', lwd = 2)
linePlotMultiEns(years = years, ens1 = fHarvest_lnd_sum_ens_wave00[without_outliers_ix_wave00,],
ens2 = fHarvest_lnd_sum_ens_wave01[without_outliers_ix_wave01,],
ens3 = fHarvest_lnd_sum_ens_wave01[level2a_ix_wave01, ],
col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
ylab = 'GtC', main = 'fHarvest')
lines(years, fHarvest_lnd_sum_stan, col = stancol, lty = 'solid', lwd = 2)
linePlotMultiEns(years = years, ens1 = treeFrac_lnd_mean_ens_wave00[without_outliers_ix_wave00,],
ens2 = treeFrac_lnd_mean_ens_wave01[without_outliers_ix_wave01,],
ens3 = treeFrac_lnd_mean_ens_wave01[level2a_ix_wave01, ],
col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
ylab = '%', main = 'treefrac'
)
lines(years, treeFrac_lnd_mean_stan, col = stancol, lty = 'solid', lwd = 2)
linePlotMultiEns(years = years, ens1 = shrubFrac_lnd_mean_ens_wave00[without_outliers_ix_wave00,],
ens2 = shrubFrac_lnd_mean_ens_wave01[without_outliers_ix_wave01,],
ens3 = shrubFrac_lnd_mean_ens_wave01[level2a_ix_wave01, ],
col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
ylab = '%', main = 'shrubfrac'
)
lines(years, shrubFrac_lnd_mean_stan, col = stancol, lty = 'solid', lwd = 2)
linePlotMultiEns(years = years, ens1 = baresoilFrac_lnd_mean_ens_wave00[without_outliers_ix_wave00,],
ens2 = baresoilFrac_lnd_mean_ens_wave01[without_outliers_ix_wave01,],
ens3 = baresoilFrac_lnd_mean_ens_wave01[level2a_ix_wave01, ],
col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
ylab = '%', main = 'baresoilfrac')
lines(years, baresoilFrac_lnd_mean_stan, col = stancol, lty = 'solid', lwd = 2)
linePlotMultiEns(years = years, c3PftFrac_lnd_mean_ens_wave00[without_outliers_ix_wave00,],
ens2 = c3PftFrac_lnd_mean_ens_wave01[without_outliers_ix_wave01,],
ens3 = c3PftFrac_lnd_mean_ens_wave01[level2a_ix_wave01, ],
col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
ylab = '%', main = 'c3PftFrac')
lines(years, c3PftFrac_lnd_mean_stan, col = stancol, lty = 'solid', lwd = 2)
linePlotMultiEns(years = years, c4PftFrac_lnd_mean_ens_wave00[without_outliers_ix_wave00,],
ens2 = c4PftFrac_lnd_mean_ens_wave01[without_outliers_ix_wave01,],
ens3 = c4PftFrac_lnd_mean_ens_wave01[level2a_ix_wave01, ],
col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
ylab = '%', main = 'c4PftFrac')
lines(years, c4PftFrac_lnd_mean_stan, col = stancol, lty = 'solid', lwd = 2)
reset()
legend('bottomright', legend = c('wave00','wave01','wave01 level2','standard'), lty = 'solid', lwd = 1.5, col = c(lcol_wave0, lcol_wave01, lcol_wave01_level2, stancol), inset = c(0.05, 0.15) )
dev.off()
null device
1
pdf(file = 'figs/fig04.pdf', width = 10, height = 12)
#pdf(file = 'figs/carbon-cycle-timeseries-anomaly-waves-constrained.pdf', width = 10, height = 12)
par(mfrow= c(3,5), las = 1, mar = c(4,4,1,0))
linePlotMultiEns(years = years, ens1 = npp_ens_anom_wave00[without_outliers_ix_wave00,],
ens2 = npp_ens_anom_wave01[without_outliers_ix_wave01,],
ens3 = npp_ens_anom_wave01[level2a_ix_wave01, ],
col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
ylab = 'GtC', main = 'NPP')
lines(years,npp_stan_anom, col = stancol, lty = 'solid', lwd = 2)
linePlotMultiEns(years = years, ens1 = nbp_ens_anom_wave00[without_outliers_ix_wave00,],
ens2 = nbp_ens_anom_wave01[without_outliers_ix_wave01,],
ens3 = nbp_ens_anom_wave01[level2a_ix_wave01, ],
col1 = lcol_wave0, col2 = lcol_wave01,col3 = lcol_wave01_level2,
ylab = 'GtC', main = 'NBP', ylim = c(-10,10))
lines(years, nbp_stan_anom, col = stancol, lty = 'solid', lwd = 2)
linePlotMultiEns(years = years, ens1 = cSoil_ens_anom_wave00[without_outliers_ix_wave00,],
ens2 = cSoil_ens_anom_wave01[without_outliers_ix_wave01,],
ens3 = cSoil_ens_anom_wave01[level2a_ix_wave01, ],
col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
ylab = 'GtC', main = 'cSoil', ylim = range(c(cSoil_ens_anom_wave00[,1], cSoil_ens_anom_wave00[,164])))
lines(years, cSoil_stan_anom, col = stancol, lty = 'solid', lwd = 2)
linePlotMultiEns(years = years, ens1 = cVeg_ens_anom_wave00[without_outliers_ix_wave00,],
ens2 = cVeg_ens_anom_wave01[without_outliers_ix_wave01,],
ens3 = cVeg_ens_anom_wave01[level2a_ix_wave01, ],
col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
ylab = 'GtC', main = 'cVeg')
lines(years, cVeg_stan_anom, col = stancol, lty = 'solid', lwd = 2)
linePlotMultiEns(years = years, ens1 = lai_lnd_mean_ens_anom_wave00[without_outliers_ix_wave00,],
ens2 = lai_lnd_mean_ens_anom_wave01[without_outliers_ix_wave01,],
ens3 = lai_lnd_mean_ens_anom_wave01[level2a_ix_wave01, ],
col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
ylab = 'GtC', main = 'Lai')
lines(years, lai_lnd_mean_stan_anom, col = stancol, lty = 'solid', lwd = 2)
linePlotMultiEns(years = years, ens1 = rh_lnd_sum_ens_anom_wave00[without_outliers_ix_wave00,],
ens2 = rh_lnd_sum_ens_anom_wave01[without_outliers_ix_wave01,],
ens3 = rh_lnd_sum_ens_anom_wave01[level2a_ix_wave01, ],
col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
ylab = 'GtC', main = 'RH')
lines(years, rh_lnd_sum_stan_anom, col = stancol, lty = 'solid', lwd = 2)
linePlotMultiEns(years = years, ens1 = fLuc_lnd_sum_ens_anom_wave00[without_outliers_ix_wave00,],
ens2 = fLuc_lnd_sum_ens_anom_wave01[without_outliers_ix_wave01,],
ens3 = fLuc_lnd_sum_ens_anom_wave01[level2a_ix_wave01, ],
col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
ylab = 'GtC', main = 'fLuc')
lines(years, fLuc_lnd_sum_stan_anom, col = stancol, lty = 'solid', lwd = 2)
linePlotMultiEns(years = years, ens1 = fHarvest_lnd_sum_ens_anom_wave00[without_outliers_ix_wave00,],
ens2 = fHarvest_lnd_sum_ens_anom_wave01[without_outliers_ix_wave01,],
ens3 = fHarvest_lnd_sum_ens_anom_wave01[level2a_ix_wave01, ],
col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
ylab = 'GtC', main = 'fHarvest')
lines(years, fHarvest_lnd_sum_stan_anom, col = stancol, lty = 'solid', lwd = 2)
linePlotMultiEns(years = years, ens1 = treeFrac_lnd_mean_ens_anom_wave00[without_outliers_ix_wave00,],
ens2 = treeFrac_lnd_mean_ens_anom_wave01[without_outliers_ix_wave01,],
ens3 = treeFrac_lnd_mean_ens_anom_wave01[level2a_ix_wave01, ],
col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
ylab = '%', main = 'treefrac'
)
lines(years, treeFrac_lnd_mean_stan_anom, col = stancol, lty = 'solid', lwd = 2)
linePlotMultiEns(years = years, ens1 = shrubFrac_lnd_mean_ens_anom_wave00[without_outliers_ix_wave00,],
ens2 = shrubFrac_lnd_mean_ens_anom_wave01[without_outliers_ix_wave01,],
ens3 = shrubFrac_lnd_mean_ens_anom_wave01[level2a_ix_wave01, ],
col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
ylab = '%', main = 'shrubfrac'
)
lines(years, shrubFrac_lnd_mean_stan_anom, col = stancol, lty = 'solid', lwd = 2)
linePlotMultiEns(years = years, ens1 = baresoilFrac_lnd_mean_ens_anom_wave00[without_outliers_ix_wave00,],
ens2 = baresoilFrac_lnd_mean_ens_anom_wave01[without_outliers_ix_wave01,],
ens3 = baresoilFrac_lnd_mean_ens_anom_wave01[level2a_ix_wave01, ],
col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
ylab = '%', main = 'baresoilfrac')
lines(years, baresoilFrac_lnd_mean_stan_anom, col = stancol, lty = 'solid', lwd = 2)
linePlotMultiEns(years = years, c3PftFrac_lnd_mean_ens_anom_wave00[without_outliers_ix_wave00,],
ens2 = c3PftFrac_lnd_mean_ens_anom_wave01[without_outliers_ix_wave01,],
ens3 = c3PftFrac_lnd_mean_ens_anom_wave01[level2a_ix_wave01, ],
col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
ylab = '%', main = 'c3PftFrac')
lines(years, c3PftFrac_lnd_mean_stan_anom, col = stancol, lty = 'solid', lwd = 2)
linePlotMultiEns(years = years, c4PftFrac_lnd_mean_ens_anom_wave00[without_outliers_ix_wave00,],
ens2 = c4PftFrac_lnd_mean_ens_anom_wave01[without_outliers_ix_wave01,],
ens3 = c4PftFrac_lnd_mean_ens_anom_wave01[level2a_ix_wave01, ],
col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
ylab = '%', main = 'c4PftFrac')
lines(years, c4PftFrac_lnd_mean_stan_anom, col = stancol, lty = 'solid', lwd = 2)
reset()
legend('bottomright', legend = c('wave00','wave01','wave01 level2','standard'), lty = 'solid', lwd = 1.5, col = c(lcol_wave0, lcol_wave01, lcol_wave01_level2, stancol), inset = c(0.05, 0.15) )
dev.off()
null device
1
Section 2.5 in Friedlingstein et al. describes how the land carbon sink is estimated.
#pdf(file = 'carbon_budget.pdf', width = 10, height = 8)
# Question: How closely should our model match this curve? Which output?
# (My guess is 'Total Land Carbon anomaly')
historical_carbon_budget <- read_excel('data/Global_Carbon_Budget_2020v1.0.xlsx', sheet = "Historical Budget", skip = 15, n_max = 270)
par(mfrow = c(3,2))
ylim = c(-1, 6)
land_sink_net <- historical_carbon_budget$`land sink` - historical_carbon_budget$`land-use change emissions`
plot(historical_carbon_budget$Year, historical_carbon_budget$`fossil emissions excluding carbonation`, main = 'fossil emissions excluding carbonation', ylab = '',
type = 'l', bty = 'n', ylim = ylim)
plot(historical_carbon_budget$Year, historical_carbon_budget$`land sink`, main = 'land sink', ylab = '', type = 'l', bty = 'n', ylim = ylim)
lines(historical_carbon_budget$Year, land_sink_net, col = 'red')
legend('topleft', legend = 'land sink - land use change emissions', col = 'red', lty = 'solid')
plot(historical_carbon_budget$Year, historical_carbon_budget$`land-use change emissions`, main = 'land use change emissions',
ylab = '', type = 'l', bty = 'n', ylim = ylim)
plot(historical_carbon_budget$Year, historical_carbon_budget$`atmospheric growth`, main = 'atmospheric growth', ylab = '',
type = 'l', bty = 'n', ylim = ylim)
plot(historical_carbon_budget$Year, historical_carbon_budget$`ocean sink`, main = 'ocean sink', type = 'l', bty = 'n', ylim = ylim)
plot(historical_carbon_budget$Year, historical_carbon_budget$`budget imbalance`, main = 'budget imbalance', ylab = '', type = 'l', bty = 'n', ylim = ylim)
#dev.off()
IPCC AR6 Chapter 5 states: “The net land carbon sink is taken as net biome productivity (NBP) and so includes any modelled net land-use change emissions. Further details on data sources and processing are available in the chapter data table (Table 5.SM.6).”
The GCP says: “The land sink is the average of several dynamic global vegetation models that reproduce the observed mean total land sink of the 1990s.”
AR6 Chapter 5 states: "The land carbon cycle components of historical ESM simulations show a larger range, with simulated cumulative land carbon uptake (1850–2014) spanning the range from –47 to +21 GtC, compared to the GCP estimate of –12 ± 50 GtC (Figure 5.23b).
#pdf(file = 'carbon_budget.pdf', width = 10, height = 8)
# Question: How closely should our model match this curve? Which output?
# (My guess is 'Total Land Carbon anomaly')
historical_carbon_budget <- read_excel('data/Global_Carbon_Budget_2020v1.0.xlsx', sheet = "Historical Budget", skip = 15, n_max = 270)
par(mfrow = c(3,2))
ylim = c(-1, 6)
land_sink_net <- historical_carbon_budget$`land sink` - historical_carbon_budget$`land-use change emissions`
plot(historical_carbon_budget$Year, historical_carbon_budget$`fossil emissions excluding carbonation`, main = 'fossil emissions excluding carbonation', ylab = '',
type = 'l', bty = 'n', ylim = ylim)
plot(historical_carbon_budget$Year, historical_carbon_budget$`land sink`, main = 'land sink', ylab = '', type = 'l', bty = 'n', ylim = ylim)
lines(historical_carbon_budget$Year, land_sink_net, col = 'red')
legend('topleft', legend = 'land sink - land use change emissions', col = 'red', lty = 'solid')
plot(historical_carbon_budget$Year, historical_carbon_budget$`land-use change emissions`, main = 'land use change emissions',
ylab = '', type = 'l', bty = 'n', ylim = ylim)
plot(historical_carbon_budget$Year, historical_carbon_budget$`atmospheric growth`, main = 'atmospheric growth', ylab = '',
type = 'l', bty = 'n', ylim = ylim)
plot(historical_carbon_budget$Year, historical_carbon_budget$`ocean sink`, main = 'ocean sink', type = 'l', bty = 'n', ylim = ylim)
plot(historical_carbon_budget$Year, historical_carbon_budget$`budget imbalance`, main = 'budget imbalance', ylab = '', type = 'l', bty = 'n', ylim = ylim)
#dev.off()
plot(historical_carbon_budget$Year, historical_carbon_budget$`land sink`,
main = 'land sink', ylab = '', type = 'l', bty = 'n', ylim = ylim)
plot(historical_carbon_budget$Year, cumsum(historical_carbon_budget$`land sink`),
main = 'land sink', ylab = '', type = 'l', bty = 'n')
#plot(historical_carbon_budget$Year, cumsum(land_sink_net, na.rm = TRUE),
# main = 'land sink', ylab = '', type = 'l', bty = 'n')
match_years_ix <- which(historical_carbon_budget$Year %in% 1850:2013)
match_years <- historical_carbon_budget$Year[match_years_ix]
cumulative_net_land_sink <- cumsum(land_sink_net[match_years_ix])
par(las = 1)
plot(match_years, cumulative_net_land_sink,
type = 'l', main = "Cumulative Net Land Sink", xlab = "", ylab = "GtC",
ylim = c(-80, 20))
par(mar = c(8,4,4,1), las = 1)
cnbp_ens_wave00 <- t(apply(nbp_ens_wave00, 1, FUN = cumsum))
cnbp_ens_wave01 <- t(apply(nbp_ens_wave01, 1, FUN = cumsum))
cnbp_stan <- cumsum(nbp_stan)
linePlotMultiEns(years = years, ens1 = cnbp_ens_wave00[without_outliers_ix_wave00,],
ens2 = cnbp_ens_wave01[without_outliers_ix_wave01,],
ens3 = cnbp_ens_wave01[level2a_ix_wave01, ],
col1 = lcol_wave0, col2 = lcol_wave01,col3 = lcol_wave01_level2,
ylab = 'Cumulative NBP (GtC)', main = 'Cumulative NBP', xlim = c(1850, 2035))
lines(years, cnbp_stan, col = 'black', lwd = 2)
lines(match_years, cumulative_net_land_sink, lty = 'dashed', col = 'black')
arrows(2022,-47, 2022,21, angle = 90, length = 0.05, code = 3)
arrows(2035,(-12-50), 2035,(-12+50), angle = 90, length = 0.05, code = 3)
text(2022, 30, 'AR6', cex = 0.8)
text(2035, 47, 'GCP', cex = 0.8)
reset()
legend('bottom', legend = c('wave00','wave01','wave01 level2','standard', 'GCP estimate'), lty = c('solid', 'solid', 'solid', 'solid', 'dashed'), lwd = 1.5, col = c(lcol_wave0, lcol_wave01, lcol_wave01_level2, stancol, stancol), horiz = FALSE, inset = 0.01, cex = 0.8)
NA
NA
selected_tags <- c('npp', 'nbp', 'cSoil', 'cVeg', 'lai_lnd_mean', 'rh_lnd_sum', 'fLuc_lnd_sum', 'fHarvest_lnd_sum', 'treeFrac_lnd_mean', 'shrubFrac_lnd_mean', 'baresoilFrac_lnd_mean', 'c3PftFrac_lnd_mean', 'c4PftFrac_lnd_mean')
# wave00
selected_tags_vec_wave00 <- paste0(selected_tags, '_ens_wave00')
mv_means_wave00 <- matrix(ncol = length(selected_tags), nrow = nrow(npp_ens_wave00))
colnames(mv_means_wave00) <- selected_tags
for(i in 1:length(selected_tags_vec_wave00)){
dat <- get(selected_tags_vec_wave00[i])
colix <- 1:ncol(dat)
trunc_dat <- dat[, tail(colix,20)]
mean_dat <- apply(trunc_dat, 1, mean, na.rm = TRUE)
mv_means_wave00[,i ] <- mean_dat
}
## wave01
selected_tags_vec_wave01 <- paste0(selected_tags, '_ens_wave01')
mv_means_wave01 <- matrix(ncol = length(selected_tags), nrow = nrow(npp_ens_wave01))
colnames(mv_means_wave01) <- selected_tags
for(i in 1:length(selected_tags_vec_wave01)){
dat <- get(selected_tags_vec_wave01[i])
colix <- 1:ncol(dat)
trunc_dat <- dat[, tail(colix,20)]
mean_dat <- apply(trunc_dat, 1, mean, na.rm = TRUE)
mv_means_wave01[,i ] <- mean_dat
}
# wave00 anomaly
selected_tags_vec_anom_wave00 <- paste0(selected_tags, '_ens_anom_wave00')
mv_means_anom_wave00 <- matrix(ncol = length(selected_tags), nrow = nrow(npp_ens_anom_wave00))
colnames(mv_means_anom_wave00) <- selected_tags
for(i in 1:length(selected_tags_vec_anom_wave00)){
dat <- get(selected_tags_vec_anom_wave00[i])
colix <- 1:ncol(dat)
trunc_dat <- dat[, tail(colix,20)]
mean_dat <- apply(trunc_dat, 1, mean, na.rm = TRUE)
mv_means_anom_wave00[,i ] <- mean_dat
}
# wave01 anomaly
selected_tags_vec_anom_wave01 <- paste0(selected_tags, '_ens_anom_wave01')
mv_means_anom_wave01 <- matrix(ncol = length(selected_tags), nrow = nrow(npp_ens_anom_wave01))
colnames(mv_means_anom_wave01) <- selected_tags
for(i in 1:length(selected_tags_vec_anom_wave01)){
dat <- get(selected_tags_vec_anom_wave01[i])
colix <- 1:ncol(dat)
trunc_dat <- dat[, tail(colix,20)]
mean_dat <- apply(trunc_dat, 1, mean, na.rm = TRUE)
mv_means_anom_wave01[,i ] <- mean_dat
}
range_wave00 <- apply(mv_means_wave00[without_outliers_ix_wave00,], 2 , range) # Wave00 sets the initial range
range_wave01 <- apply(mv_means_wave01[without_outliers_ix_wave01,], 2, range)
range_wave01_level2a <- apply(mv_means_wave01[level2a_ix_wave01, ], 2, range)
# What is the output range of the level2a modern value, expressed as a proportion of the initial ensemble?
range_prop_wave01_level2a <- (apply(range_wave01_level2a, 2, diff) / apply(range_wave00, 2, diff)) *100
range_anom_wave00 <- apply(mv_means_anom_wave00[without_outliers_ix_wave00,], 2 , range) # Wave00 sets the initial range
range_anom_wave01 <- apply(mv_means_anom_wave01[without_outliers_ix_wave01,], 2, range)
range_anom_wave01_level2a <- apply(mv_means_anom_wave01[level2a_ix_wave01, ], 2, range)
# What is the output range of the level2a modern value, expressed as a proportion of the initial ensemble?
range_prop_anom_wave01_level2a <- (apply(range_anom_wave01_level2a, 2, diff) / apply(range_anom_wave00, 2, diff)) *100
Proportion of the initial range of model output that is covered by level 2 constrained ensemble (%):
all_waves_mv <- rbind(mv_means_wave00[without_outliers_ix_wave00,], mv_means_wave01[without_outliers_ix_wave01,])
range_all_waves <- apply(all_waves_mv, 2, range)
level2_ix_all_waves_mv <- which(all_waves_mv[,'nbp'] > 0 &
all_waves_mv[,'npp'] > 35 &
all_waves_mv[,'npp'] < 80 &
all_waves_mv[,'cSoil'] > 750 &
all_waves_mv[,'cSoil'] < 3000 &
all_waves_mv[,'cVeg'] > 300 &
all_waves_mv[,'cVeg'] < 800
)
all_waves_mv_level2 <- all_waves_mv[level2_ix_all_waves_mv, ]
all_waves_col <- rep(zred, nrow(all_waves_mv))
all_waves_level2_col <- all_waves_col
all_waves_level2_col[level2_ix_all_waves_mv] <- zblue
range_all_waves_level2 <- apply(all_waves_mv[level2_ix_all_waves_mv, ], 2, range)
range_prop_all_waves_level2 <-(apply(range_all_waves_level2, 2, diff) / apply(range_all_waves, 2, diff)) *100
# Anomaly
all_waves_anom_mv <- rbind(mv_means_anom_wave00[without_outliers_ix_wave00,], mv_means_anom_wave01[without_outliers_ix_wave01,])
range_all_waves_anom <- apply(all_waves_anom_mv, 2, range)
all_waves_anom_mv_level2 <- all_waves_anom_mv[level2_ix_all_waves_mv, ]
range_all_waves_anom_level2 <- apply(all_waves_anom_mv[level2_ix_all_waves_mv, ], 2, range)
range_prop_all_waves_anom_level2 <-(apply(range_all_waves_anom_level2, 2, diff) / apply(range_all_waves_anom, 2, diff)) *100
print(round(range_prop_all_waves_level2,1))
npp nbp cSoil cVeg lai_lnd_mean rh_lnd_sum
26.6 55.1 48.3 23.5 60.9 26.0
fLuc_lnd_sum fHarvest_lnd_sum treeFrac_lnd_mean shrubFrac_lnd_mean baresoilFrac_lnd_mean c3PftFrac_lnd_mean
51.5 60.7 88.3 85.4 52.9 58.0
c4PftFrac_lnd_mean
43.9
pdf(file = "figs/induced_constraint_barplot.pdf")
range_prop_round = round(range_prop_all_waves_level2,1)
range_prop_round_mat <- rbind(range_prop_round, (100 - range_prop_round))
range_prop_anom_round = round(range_prop_all_waves_anom_level2,1)
range_prop_anom_round_mat <- rbind(range_prop_anom_round, (100 - range_prop_anom_round))
par(las = 1, mar = c(8,1,3,2), mfrow = c(1,2), oma = c(1,10, 1,1))
b1 <- barplot(t(apply(range_prop_round_mat,1,rev)),
horiz = TRUE,
col = c( zblue, zred),
xlab = 'Ensemble output range (%)',
main = 'Absolute'
)
text(0, rev(b1), range_prop_round, pos = 4, cex = 0.9)
b2 <- barplot(t(apply(range_prop_anom_round_mat,1,rev)),
horiz = TRUE,
col = c( zblue, zred),
xlab = 'Ensemble output range (%)',
names.arg = rep('', 13),
main = 'Anomaly'
)
text(0, rev(b2), range_prop_anom_round, pos = 4, cex = 0.9)
axis(1)
reset()
legend('bottom', legend = c('NROY', 'Ruled out'), fill = c(zblue, zred), horiz = TRUE, inset = 0.02 )
dev.off()
null device
1
axat = length(range_prop_round):1
par(las = 1, mar = c(5,10,3,2), xaxs = 'i')
plot(rev(range_prop_round), axat, xlim = c(0,100), axes = FALSE, type = 'n', ylab = '', xlab = 'Level 2 NROY proportion of whole range (%)')
#abline(h = axat, col = 'grey', lty = 'dashed')
#abline(v = 100)
points(rev(range_prop_anom_round), rev(axat), col = 'grey', pch = 19)
segments(x0 = 0, y0 = rev(axat), x1 = rev(range_prop_anom_round), y1 = rev(axat), col = 'grey', lty = 'dashed')
points(rev(range_prop_round), rev(axat), col = 'black', pch = 19)
segments(x0 = 0, y0 = rev(axat), x1 = rev(range_prop_round), y1 = rev(axat))
#segments(x0 = 0, y0 = rev(axat), x1 = rev(range_prop_anom_round), y1 = rev(axat), col = 'grey', lty = 'dashed')
axis(2, at = axat, labels = selected_tags)
axis(1)
reset()
legend('top', c('Absolute', 'Anomaly'), pch = 19, col = c('black', 'grey'), horiz = TRUE, inset = 0.03)
NA
NA
#pdf(file = 'figs/output_pairs.pdf', width = 10, height = 10)
pairs(all_waves_mv, col = all_waves_level2_col, lower.panel = NULL, gap = 0,
pch = 20,
labels = 1:13
)
reset()
legend('topleft', legend = c('NROY', 'Ruled out' ), pch = 20, col = c(zblue, zred), inset = c(0.03, 0.25))
legend('left', legend = c(paste(1:13, colnames(all_waves_mv))), inset = 0.03, cex = 1)
#dev.off()
# Can this go in common data? Would be needed for checking emulator fits
# Create fit lists for the combined data wave00 level 1a and wave01
Y_const_level1a_wave01_scaled_list <- mat2list(Y_const_level1a_wave01_scaled)
fit_list_const_level1a_wave01 <- mclapply(X = Y_const_level1a_wave01_scaled_list , FUN = km, formula = ~., design = X_level1a_wave01,
mc.cores = 4, control = list(trace = FALSE))
Y_const_level1a_wave01_scaled_pred <- multiPred(Y = Y_const_level1a_wave01_scaled, Xpred = X_unif, fit_list = fit_list_const_level1a_wave01)
level2_ix_em_unif_wave00_wave01 <- which(Y_const_level1a_wave01_scaled_pred$pred_mean[,'nbp_lnd_sum'] > 0 &
Y_const_level1a_wave01_scaled_pred$pred_mean[,'npp_nlim_lnd_sum'] > 35 &
Y_const_level1a_wave01_scaled_pred$pred_mean[,'npp_nlim_lnd_sum'] < 80 &
Y_const_level1a_wave01_scaled_pred$pred_mean[,'cSoil_lnd_sum'] > 750 &
Y_const_level1a_wave01_scaled_pred$pred_mean[,'cSoil_lnd_sum'] < 3000 &
Y_const_level1a_wave01_scaled_pred$pred_mean[,'cVeg_lnd_sum'] > 300 &
Y_const_level1a_wave01_scaled_pred$pred_mean[,'cVeg_lnd_sum'] < 800
)
(length(level2_ix_em_unif_wave00_wave01) / nunif) * 100
X_stan_norm <- normalize(matrix(rep(1, 32), nrow = 1), wrt = lhs)
colnames(X_unif) <- 1:32
pdf(file = 'figs/fig08.pdf', width = 12, height = 12)
#pdf(file = 'figs/pairs_level2_ix_em_unif_wave00_wave01.pdf', width = 12, height = 12)
par(oma = c(0,0,0,3), bg = 'white')
panel_hist_local <- function(x, ...)
{
usr <- par("usr"); on.exit(par(usr))
par(usr = c(usr[1:2], 0, 1.5) )
h <- hist(x, plot = FALSE)
breaks <- h$breaks; nB <- length(breaks)
y <- h$counts; y <- y/max(y)
rect(breaks[-nB], 0, breaks[-1], y, col = "cyan", ...)
}
pairs(rbind(X_unif[level2_ix_em_unif_wave00_wave01, ], X_stan_norm),
gap = 0, lower.panel = NULL, xlim = c(0,1), ylim = c(0,1),
panel = dfunc_up_truth,
diag.panel = panel_hist_local,
cex.labels = 1,
col.axis = 'white',
dfunc_col = rb
)
image.plot(legend.only = TRUE,
zlim = c(0,1),
col = rb,
legend.args = list(text = 'Density of model runs matching the criteria', side = 3, line = 1),
legend.shrink = 0.6,
horizontal = TRUE
)
legend('left', legend = paste(1:32, colnames(lhs)), cex = 1, bty = 'n')
dev.off()